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Abstract 


A new probabilistic facility-centric approach to lightning strike location has been developed. 

This process uses the bivariate Gaussian distribution of probability density provided by the current 
lightning location error ellipse for the most likely location of a lightning stroke and integrates it to 
determine the probability that the stroke is inside any specified radius of any location, even if that location 
is not centered on or even with the location error ellipse. This technique is adapted from a method of 
calculating the probability of debris collision .with spacecraft. Such a technique is important in spaceport 
processing activities because it allows engineers to quantify the risk of induced current damage to critical 
electronics due to nearby lightning strokes. This technique was tested extensively and is now in use by 
space launch organizations at Kennedy Space Center and Cape Canaveral Air Force Station. Future 
applications could include forensic meteorology. 
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Nomenclature 


0 = angle between the collision plane coordinates and the rotated collision plane coordinates 

Pxz = correlation coefficient of x and z 

o K ' = standard deviation of x-coordinate of the diagonalized covariance ellipse in the rotated 
coordinate system 

o H ’ — standard deviation of z-coordinate of the diagonalized covariance ellipse in the rotated 
coordinate system 
g x - standard deviation of x 

o z — standard deviation of z 

o x > = standard deviation of x in the rotated coordinate system 

a z ’ - standard deviation of z in the rotated coordinate system 

ji K = x-coordinate of target circle in the (X’, Z’) coordinate system 

ji H - z-coordinate of target circle in the (X’, Z’) coordinate system 

A = collision cross-sectional area (nautical miles 2 , nmi 2 ) 

d0 = the angle between two points on the target ellipse 

dH = integration step 

H = intermediate variable in the lightning probability algorithm 

P - probability 

pdf = probability distribution function 

r A = radius of circle of cross-sectional area, A (nautical miles, nmi) 

R = radius of ellipse (nmi) 

R1 = distance to first point on target ellipse (nmi) 

R2 = distance to second point on target ellipse (nmi) 

x = horizontal rectangular coordinate in the collision plane (nautical miles, nmi) 

x 5 = horizontal rectangular coordinate in the rotated collision plane (nmi) 

x 55 = transformation variable that circularizes the x-component of the probability ellipse (nmi) 

x e = nominal distance of closed approach of two colliding objects (nmi) 

W = intermediate variable in the lightning probability algorithm 

z = vertical rectangular coordinate in the collision plane (nautical miles, nmi) 

z ' = vertical rectangular coordinate in the collision plane (nautical miles, nmi) 

z” = transformation variable that circularizes the z-component of the probability ellipse (nmi) 
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Introduction 


The ability to a Introduction accurately estimate the probability that an individual nearby cloud- 
to-ground lightning stroke was within a specified distance of any specified spaceport processing facility at 
Kennedy Space Center (KSC) or Cape Canaveral Air Force Station (CCAFS) is important to processing 
payloads and space launch vehicles before launch. Such estimates allow engineers to decide if inspection 
of electronics systems aboard satellite payloads, space launch vehicles, and ground support equipment is 
warranted due to induced currents from that stroke. If induced current damage has occurred, inspections 
of the electronics are critical to identify required fixes and avoid degraded performance or failure of the 
satellite or space launch vehicle. However, inspections are costly both financially and in terms of delayed 
processing for space launch activities. As such, it is important these inspections be avoided if not needed. 
At KSC/CCAFS, one of the main purposes of the Four Dimensional Lightning Surveillance System 
(4DLSS) (Murphy, 2008, Roeder, 2010) is detection of nearby strokes and determination of their peak 
current to support decisions to inspect electronics (Flinn, 2010a, Flinn, 2010b, Roeder, 2005). The high 
frequency of lightning occurrence in East Central Florida combined with the large amount of complex 
sensitive electronics in satellite payloads, space launch vehicles, and associated facilities makes those 
decisions critically important to space launch processing. The 4DLSS provides the data for 50th 
percentile location error ellipses for the best location for each stroke, which is then scaled to 95th or 99th 
percentile ellipses depending on customer requirements. This error ellipse is necessarily centered on the 
best location of the lightning stroke. The 4DLSS, however, has not been able to provide the probability 
of the stroke being within a customer specified distance of a point of interest. This paper presents a new 
method to convert the 4DLSS 50th percentile location error ellipse for best location of any stroke into the 
probability that the stroke was within any radius of any facility at CCAFS/KSC. This technique could be 
adapted for use with National Lightning Detection Network (NLDN) data. This new probabilistic facility- 
centric technique is a significant improvement over the stroke-centric location error ellipses the 45th 
Weather Squadron (45 WS) has provided in the past. This technique is adapted from a method of 
calculating the probability of debris collision with spacecraft (Chan, 2008, Leleux, 2002, Patera, 2001). 


Background 


Methodology 


In spacecraft collision probability and other applications, at the instant of “nominal” closest 
approach, the position uncertainty of the collision object relative to the asset being protected is described 
by a bivariate Gaussian probability density function (pdf) (Chan, 2008, Patera, 2001, Alfano, 2006, 
Alfano 2007), as shown in the following equation 


/(*»*) = 


2 m x <r^l-p 2 x 


f f- r- M 4- 

<T t <7 


/ 2(1— 


( 1 ) 


where o x and c z = the standard deviations of x and z, p xz = correlation coefficient of x and z, x and z are 
the designations for the rectangular coordinates in the collision plane. 

The probability of collision (Eq. 2) is given by the two-dimensional integral, where A is the 
collision cross-sectional area which is a circle with radius, r A (Chan, 2008). 
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( 2 ) 


r= JJ/fe z)dxdz 


A 


There is no known analytical solution to the above integral when the two standard deviations o x 
and o z are not equal. The solution is found by performing a numerical integration of the two dimensional 
Gaussian pdf (Chan, 2008, Patera, 2001, Alfano, 2006, Alfano 2007). 

The geometry used for spaceflight collision probability can also be used for estimation of the 
probability of an individual nearby lightning stroke contacting the surface within a specified distance of a 
specified point of interest as shown in Fig. 1. In Fig. 1, a is the heading of the semi-major axis of the 
lightning location uncertainty ellipse from true north and 0 is the angle between the semi-major axis of 
the lightning location uncertainty ellipse and line connecting the center of the lightning uncertainty ellipse 
and the center of the area of interest. Two methods of integrating the above probability were tested and 
gave identical results. The first solution method was based on an algorithm by Patera (2001) and the 
second solution method was based on an algorithm by Chan (201 1). Chan’s algorithm ran much faster 
and therefore was selected as the algorithm for the 45 WS lightning probability program. 


The Numerical Integration Technique (Chan, 2011) 

This numerical integration technique is one in which the miss distance is given by a non-central 
chi distribution with unequal variances (Chan, 201 1). The covariance matrix corresponding to the 
bivariate Gaussian pdf in Eq. 1 is (Chan, 2008): 


Area around point of interest 


l-o uncertain! 
around lightnii 



imple 
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(3) 


Pxz°x°z 

When the correlation coefficient, p^, is not zero, there are undesirable off-diagonal terms that 
overly complicate the calculation. In order to eliminate these terms, the coordinate system (x,z) is rotated 
to a new coordinate system (x’,z’) such that the major and minor axes of the ellipse associated with the 
covariance are aligned along the coordinate axes and the new covariance matrix is (Chan, 2008): 






The angle, 0, between the two coordinate systems is (Chan, 2008): 


(4) 


d = — tan 1 
2 


2 P M <* x <r t 


(5) 


The KSC/CCAFS 4DLSS system does not provide the covariance matrix, but instead provides 
the semi-major axis, semi-minor axis, and the orientation of the semi-major axis of the 50% location error 
ellipse relative to north. Therefore the angle, 0, in Eq. 5 is found using geometry where 0 is the angle 
between the semi-major axis of the lightning location uncertainty ellipse and line connecting the center of 
the lightning uncertainty ellipse and the center of the area of interest. 

In the (x J , z’) system, the Eq. 1 pdf becomes (Chan, 2008) 





2/T<7 x ,< 7 z , 


and Eq. 2, the collision probability becomes (Chan, 2008) 


( 6 ) 


P = - 


1 


2na,o. 



J dx'dz ' 


(7) 


where 


A’= A,r ., = r A ,x'„ = x cos 6,z'-x sin# (8) 

A A p e p e y ’ 

For spacecraft collision, x e is the nominal distance of closest approach of the two colliding objects 
and (x’ p , z’ p ) are the coordinates of the spacecraft relative to the debris. For lightning strike probability, 
x e (the distance between the position of the center of the strike location ellipse and the position of the 
target area) is calculated using the Haversine distance formula. 
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The standard deviations in the new rotated coordinate system are calculated by dividing the semi- 
major and semi -minor axes of the 50% lightning positional confidence ellipse by the scaling constant 
used to scale standard error to the 50% confidence level. The scaling constant is: 

k = V-2*ln(l-0.50) (9) 


The probability is given by (Chan, 2011) 


VFr 




ijlncr 


-(H-Mh fuel + e -(H+VH file] 


1 erf(Z t )+erf(Z 2 )]lH 


H 0 


where 


2<r t 


Z 2 = y(^-// 2 ) + //J/V2a 


( 10 ) 


( 11 ) 


The parameters ^ K and ^ H are the coordinates of the target circle in the (X’, Z’) coordinate 
system; and ® K and ® H are the standard deviations of the diagonalized covariance ellipse shown in Eq. 
4. The derivation of equations (10) and (11) above is shown in further detail in (Chan, 2011). A detailed 
example of the calculations using a real-world case is provided in Appendix-A. The Excel Visual Basic 
code to implement these calculations is shown in Appendix-B. 


Evaluation 

The probability that any lightning strike is within any radius of any point of interest would be 
extremely difficult to estimate intuitively. As a result, given the high impact of the decisions on space 
launch operations, the tool developed for this application was extensively tested. Tests were conducted 
and are discussed in the following sections: 1) known mathematical solutions, and 2) examination of 
real-world events. The new technique passed all of the tests. Tests were also conducted to assure the 
probabilities calculated using the algorithm of Chan (201 1) matched probabilities calculated using the 
algorithm of Patera (2001). The probabilities calculated from the two algorithms were identical and thus 
are not shown here. 


Test Set 1 

The first set of testing compared the lightning strike probability calculated using the 45 WS 
lightning strike spreadsheet (which uses an adaptation of the numerical integration algorithm by Chan, 

201 1, to the corresponding circular probability from the CRC Handbook of Tables for Probability and 
Statistics (Beyer, 1968). Table 1 shows the probability from the new numerical integration technique for 
various inputs and the corresponding correct probability from the CRC Handbook. The values matched to 
within a tenth of a percent. These errors in the final digit may be due to round-off error. 
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Table 1. Calculated probability vs. CRC Handbook probability for various inputs 


Semi- 

major 

axis 

(nmi) 

Semi- 

minor 

axis 

(nmi) 

Heading 
of semi- 
major 
axis from 
true 
North 

Point 

Of 

Interest 

latitude 

Point 

Of 

Interest 

long- 

itude 

Strike 

Latitude 

Strike 

Long- 

itude 

Radius 
around 
Point Of 
Interest 
(nmi) 

Calcu- 

lated 

prob- 

ability 

CRC 

Hand- 

book 

prob- 

ability 

(Beyer, 

1968) 

3 

3 

15 

28.6082 

-80.6041 

28.6995 

-80.6041 

3 

0.095 

0.095 

3 

3 

15 

28.6082 

-80.6041 

28.631 

-80.6041 

3 

0.453 

0.452 

3 

3 

15 

28.6082 

-80.6041 

28.608 

-80.6041 

3 

0.500 

0.499 

1 

1 

15 

28.6082 

-80.6041 

28.608 

-80.6041 

1 

0.500 

0.499 

1 

1 

15 

28.6082 

-80.6041 

28.631 

-80.6041 

1 

0.200 

0.200 

1 

1 

15 

28.6082 

-80.6041 

28.6995 

-80.6041 

1 

0.000 

0.000 

1 

1 

15 

28.6082 

-80.6041 

28.608 

-80.6041 

2 

0.937 

0.938 


Test Set 2 

The second type of testing analyzed six real-world lightning strikes near Space Launch Complex 
39A on 3 August 2009. Figure 2 shows the spreadsheet used to generate the lightning report for those six 
strikes. Additional data on three of these six strikes are in Table 2. These strikes were selected because the 
closest point on the lightning position uncertainty ellipse was within 0.45 nautical miles of Launch 
Complex 39A, which is the key radius for assessing the need to inspect electronics for induced current 
damage to the Space Shuttle. Figures 3 through 5 are Google Maps depictions of three of these six 
strokes. In Figures 3-5, the black ellipse is the 99% lightning location uncertainty ellipse. The white 
circle is a 0.45 nmi radius around the point of interest. The probabilities for a small area around a facility, 
even for a nearby stroke, may appear to be surprisingly low. For example, figures 3 and 4 respectively 
illustrate a 53.8 % and 7.7% probability that the lightning strike occurred within the area of interest 
while figure 5 shows that a strike just 0.65 nautical miles away had only a 1.1% probability of being 
within the 0.45 nautical mile radius of Launch Complex 39A. All calculated probabilities are consistent 
with these real-world events. 
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Figure 2. Sample of lightning strikes where the closest point on the lightning position uncertainty 
ellipse was within 0.45 nmi of Launch Complex 39A on 3 August 2009. 
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The KSC Electromagnetic Environmental Effects (EEE) Panel requested six more real-world 
lightning strikes be investigated. These were recently investigated lightning strikes near Launch 
Complexes 39A or 39B where there was camera verification of the location of the strike. The EEE Panel 
wanted to compare the results of the new facility-centric probabilistic technique to these cases where the 
true answers were known unambiguously. The data used for this analysis are in Table 3. Both 4DLSS 
and National Lightning Data Network (NLDN) cases were examined, depending upon which sensor 
system recorded the stroke. CGLSS strokes were obtained from 45WS 4DLSS. The NLDN usually 
provided flash data, so NLDN return stroke data were purchased as special StrikeNet 1 reports from 
Vaisala Corporation. This was done to match the return strokes routinely provided by 4DLSS. Figures 6 
through 8 show the probability results from these cases. In Figures 6-8, the black ellipse is the 99% 
lightning location uncertainty ellipse. The white circle is a 0.45 nmi radius around the point of interest. 
As with the previous real-world tests, all calculated probabilities were consistent with these additional 
real-world events. 

Table 2. Input values used for scenarios shown in Figures 3 through 5. 


Figure 

Semi- 

Semi- 

Confidence 

Heading 

Point of 

Point of 

Strike 

Strike 

Radius 


major axis 

major axis 


(from 

interest 

interest 

latitude 

longitude 

around 


of 50% 

of 50% 


true 

latitude 

longitude 

(°N) 

CW) 

point 


confidence 

confidence 


North) 

(°N) 

(°W) 



of 


ellipse 

ellipse 


of semi- 





interest 


(km) 

(km) 


major 





(nmi) 





axis 






3 

0.4 

0.2 

0.99 

300.7 

28.60827 

-80.6041 

28.6114 

-80.6113 

0.45 

4 

0.3 

0.2 

0.99 

293 

28.60827 

-80.6041 

28.6178 

-80.6069 

0.45 

5 

0.2 

0.1 

0.99 

20.3 

28.60827 

-80.6041 

28.5995 

-80.6113 

0.45 


1 Vaisala, Inc., 2006, “Vaisala StrikeNet Information,” URL: http://www.vaisala.com/files/StrikeNet-Brochure.pdf 
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Figure 3. Google Maps visualization of the 99% 
confidence uncertainty ellipse for one of the closest 
lightning strikes to Complex 39A on 03 August 
2009. 



Figure 4. Google Maps visualization of the 99% 
confidence uncertainty ellipse for a lightning strike 
near Complex 39 A on 03 August 2009. 
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Figure 5. Google Maps visualization of the 99% 
confidence uncertainty ellipse for nearby lightning 
strike to Complex 39A on 03 August 2009. 


Table 3. Input values used for scenarios shown in Figures 6 through 8. 


Figure 

Semi- 
major axis 
of 50% 
confidence 
ellipse 
(km) 

Semi- 
major axis 
of 50% 
confidence 
ellipse 
(km) 

Confidence 

Heading 
(from 
true 
North) 
of semi- 
major 
axis 

Point of 
interest 
latitude 
(°N) 

Point of 
interest 
longitude 
(°W) 

Strike 

latitude 

(°N) 

Strike 

longitude 

(°W) 

Radius 

around 

point 

of 

interest 

(nmi) 

6 

0.6 

0.4 

0.99 

82 

28.60827 

-80.6041 

28.6069 

-80.6087 

0.45 

7 

0.4 

0.4 

0.99 

95 

28.60827 

-80.6041 

28.6057 

-80.6085 

0.45 

8 

0.2 

0.1 

0.99 

72 

28.62716 

-80.6275 

28.6275 

-80.6202 

0.45 
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Figure 6. Illustrates a probability of 69.1% of a 
lightning strike of amplitude -43.0 kA detected by 
NLDN occurring 0.26 nmi from the center of 
Launch Complex 39A on 8/16/2009. 



Figure 7. Illustrates a probability of 74.7% of a 
lightning strike of amplitude -71 .4 kA detected by 
NLDN occurring 0.28 nautical miles from the 
center of Launch Complex 39 A on 10/14/2009. 
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Figure 8. Illustrates a probability of 99.9996% of a 
lightning strike of amplitude -21 .7 kA detected by 
CGLSS occurring 0.04 nmi from the center of 
Launch Complex 39B on 6/27/2009. 


Other Applications 


The techniques and methods described in this paper clearly have application reaching far beyond 
the space program uses for which they were designed. The list of potential applications is many and 
varied and would be of interest to anyone seeking information pertaining to probability of lightning strike 
locations, such as the power industry, aviation, or any industry sensitive to electrical overloads. This 
methodology is also applicable to forensic meteorology [13], where the question of whether lightning 
struck at or near a particular location is an issue in litigation. This technique could be applied to NLDN 
or any other cloud-to-ground lightning detection system that provides location error ellipses. 


Conclusion 

A probabilistic facility-centric approach to lightning location has been developed to calculate the 
probability that any nearby cloud-to-ground lightning stroke occurred within any radius of any point of 
interest. In practice, this provides the probability that a nearby lightning stroke was within a key distance 
of a facility, rather than within the error ellipses centered on the stroke. This process uses the bivariate 
Gaussian distribution of probability density provided by the current lightning location error ellipse for the 
most likely location of a lightning stroke and integrates it to determine the probability that the stroke is 
inside any specified radius. This new facility-centric technique was tested extensively and is much more 
useful to the space launch customers. 
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Appendix A - Example Lightning Probability Calculation 

This appendix is an example of calculating the probability of any lightning stroke with a known 
error ellipse being within a circle of any radius around any point. It is provided to clarify the calculation 
process. An example calculation is shown in Table 4. 

This example is a real-world event from a lightning strike near .the Space Shuttle launch pad 39A 
at 02:35 GMT on 16 August 2009 (ref. Figure 6). Although the lightning data usually are from the cloud- 
to-ground component of the Four Dimensional Lightning Surveillance System (CG-4DLSS) (Murphy et 
al 2008, Roeder 2010) in this example a lightning stroke from the NLDN is used. We sometimes use 
StrikeNet reports that provide stroke data from the NLDN to double check the CG-4DLSS report. The 
lightning stroke input values are shown in row 1 of Table 3. 

The following assumptions are applicable to the example calculation in Table 4. The location of 
Launch Pad 39A is 28.60827486 north latitude (or 0.499309 radians) and 80.6041 1653 west longitude (or 
-1 .406807 radians). This is also the center of the circle in which the lightning probability will be 
calculated. The desired radius for probability of lightning around Launch Pad 39A is 0.45 nautical miles 
(nmi). This lightning stroke occurred on 16 August 2009 at 0235 GMT. 

Table 4. Lightning strike probability calculation process. 

Step Action Equation and other information Calculation and Result 


Convert semi- 

1 nmi = 1.852 km 

0.6 km = 0.324 nmi semi- 

major and 


major axis 

semi-minor 


0.4 km = 0.216 nmi semi- 

axes from km 


minor axis 

to nmi 
Calculate 

Haversine Distance Formula: 

Distance = 3443.920086 * 

distance from 

Distance = Earth Radius * C 

7.42 17x1 0‘ 5 

lightning stroke 

• Earth Radius = 3443.920086 nmi 

' — 0.2556 nmi 

to center of 

• C = 2*Atn2[Jl-A,JJ) 


circle 

Atn2 is a two parameter arc tangent 

C = 2*Atn2{sqr(l - 
1.377x1 O’ 9 ), 


function which returns values in all four 


quadrants. 

sqr(1.377xl0’ 9 )} 
= 7.4217xl0' 5 


• A = sin(dlat/2)*sin(dlat/2) + cos( 


(target lat))* cos((stroke lat)) 

A = sin(l. 200x1 0' 5 /2)* 


*sin(dlon/2)*sin(dlon/2) 


• dlat = latitude difference from 

sin(l. 200 x1 0' 5 /2)* 


target to stroke = 

+ 


28.60827° - 28.6069° = 

cos(0.4993)*cos(0.4993)* 


2.39959xl0' 5 (radians) 

sin(4.000xl0' 5 /2) 


• dlon = longitude difference from 

*sin(4. 000x1 0‘ 5 /2) 


circle to stroke = -80.6041 1 °- 

= 1.3770x1 0' 9 

Perform 

80.6085° = 7.99967X 10‘ 5 
(radians) 

• X = (Longitude of Target - Longitude 

X = (-1.4068- (-1.4069) * 

coordinate 

of Stroke ) * Cos (Latitude of Strike) 

Cos (0.4993) 


system rotation 
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to eliminate the 
off-diagonal 
term in the 
covariance 
matrix. 


4 Calculate the 
standard 
deviations in 
the new rotated 
coordinate 
system. 


5 Calculate the 
probability that 
lightning stroke 
was within the 
target area of 
interest by 
performing a 
numerical 
integration 
using 

Simpson’s rule 
of the lightning 
uncertainty 
ellipse over the 
area of the 
circle around 
the target of 
interest. 


• Z = Latitude of Target - Latitude of 
Stroke 

• 0 - a - ((tc/ 2) - Atn2(X,Z) 

• a is the orientation angle of the 50% 
lightning positional confidence ellipse 

• X’ = miss distance*Cos(0 (coordinate 
system rotation angle)) 

• Z’ = miss distance*Sin(0) 


• Gx> = Semi-major axis of the 50% 
lightning positional confidence ellipse 
/elliptical scaling constant used to scale 
standard error to the 50% confidence 
level 

• o z > = Semi-minor axis of the 50% 
lightning positional confidence ellipse 
/elliptical scaling constant used to scale 
standard error to the 50% confidence 
level 

• Elliptical scaling constant, k, is: 

2 * ln(l — probability) 


2y[27T< 


1 4W\ 

- r 

r H 0 L 


Zi = 


e [ H Vh) /2o h [ h+/i h) /2o h 


rfiz^+erfizJjH 

yKW-H^+flA/JlCTj, 


• p K and li 1 1 are the coordinates of the target 
circle in the (X’, Z’) coordinate system 

• o K and o H are equal to Gx> and O/- which are 
the standard deviations of the diagonalized 
covariance matrix. 

• W = Radius around target 2 

• N = the number of iterations to perform 
in the integration (for this example, N is 
set to 200). 

• DH = -Jw / N = integration step 

• B, C, and D are intermediate variables in 


= 7.0231 XI 0' 5 ” 

Z = 0.49931 -0.49929 
= 2.400 X 10' 5 

0 = 1.431 - 1.571 - 
Atn2(7.023 X 10 -5 , 2.400 
X 10' 5 ) =0.1896 

X’ = 0.2556*Cos(0.1896) 
= 0.2510 

Z’= 0.2556*Sin(0.1896) 
= 0.0482 

k = 2 * ln(i — 0.50) = 

1.1774 

o x . = 0.3240/1.1774 = 
0.2752 

o z . = 0.2160/1.1774 = 
0.1834 


W = Radius around target 2 
W = 0.45 2 = 0.2025 

DH = 4w / N 
DH = V0.2025 /200 = 
0.00225 

B = 42 * o x » 

B= 4l *0.2752 = 0.3891 

C = X’/ B 

C = 0.2510/0.3891 = 
0.6451 

D = 1 / (2 * yfljr * o z >) 

D = 1/(2 * Jin* 0.1834) = 
1.0874 

H = iteration no. * DH 
H= 199*0.00225 
H = 0.4478 
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the algorithm corresponding to various 
parts of the probability equation shown 
above 

• B = V2 * cx’ 

• C = X7B 

• D = 1/(2 * 42k * o z .) 

• A, H, zl, z2, E, F, Erf(zl), Erf(z2), Q and 
sum are intermediate variables in the 
algorithm corresponding to various parts 
of the probability equation shown above 

• A loop is performed for j = 1 to 199. 

This example is shown for j = 199. 



zl =A/B-C 

zl =4.494 X 10' 2 / 0.3891 

-0.6451 

zl = -0.5296 

z2 = A / B + C 

Z 2 = 4.494 X 10' 2 / 0.3891 

+ 0.6451 

z2 = 0.7606 


• Sum = 0 

• Begin Loop here: H = j* DH 

• A= 

• zl =A/B-C 

• z2 = A/B + C 

• Erf(x) = error function = 

—== \e dt 
y 7T 0 

• E = (H - Z’) 2 / (2 * oz 2 ) 

• F = (H + Z’) 2 / (2 * o z 2 ) 

• Q(j) = (e' E + e' F ) * (Erf(Zl) + 
Erf(Z2)) 

• sum = sum + (3 - (-1) j ) * Q(j) 


Erf(Zl) = 
ErrorFunction(z 1 ) 
Erf(-0.5296) = -0.5461 

Erf(Z2) = 
ErrorFunction(z2) 
Erf(0.7606) = 0.7179 
E = (H - Z’) 2 / (2 * o z 2 ) 
E = (0.4478 - 
0.0482) 2 /(2*0. 1 834 2 ) 

E = 2.372 

F = (H + Z’) 2 / (2 * o z 2 ) 
F = (0.4478 + 
0.0482) 2 /(2*0.1834 2 ) 

F = 3.654 


• sum = sum + Q(0) + Q(N) QG ) = (e‘ E + e F ) * 

• Probability = D*sum*DH/3 (Erf(Zl) + Erf(Z2)) 

Q(199) = (e 2372 + e' 3 ' 654 )* 
(-0.5461+0.7179) = 
2.0467 X 1 O’ 2 


sum = sum + (3 - (-1) J ) * 

QG) 

sum = sum + (3 - (-1) '") 
*2.047 XI O' 2 
sum = 844.8952 
End Loop 


sum = sum + Q(0) + Q(N) 
sum = 844.8952 + 2.9361 
+ 0 = 847.8317 


Probability = D * sum * 
DH/3 

Probability = 1.0874 * 
847.8317 *0.00225/3 = 
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0.6914 


E ="(H - Z’) 2 / (2 * o z 2 ) 

E = (0.4478 - 

0.04 82) 2 /(2*0. 1 834 2 ) 

E = 2.372 

F = (H + Z’) 2 /(2 *g z ’ 2 ) 

F = (0.4478 + 
0.0482) 2 /(2*0. 1 834 2 ) 

F = 3.654 

QG) = (e E + e- F )* 
(Erf(Zl) + Erf(Z2)) 
Q(199) = (e 2372 + e' 3 ' 654 )* 
(-0.5461+0.7179) = 
2.0467 X 1 0‘ 2 

sum = sum + (3 - (-1) J ) * 

QG) 

sum = sum + (3 - (-1) 199 ) 
* 2.047 X 10‘ 2 
sum = 844.8952 
End Loop 


sum = sum + Q(0) + Q(N) 
sum = 844.8952 + 2.9361 
+ 0 = 847.8317 


Probability = D * sum * 
DH/3 

Probability = 1 .0874 * 
847.8317 *0.00225/3 = 
0.6914 
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Appendix B - Excel VBA Macro for 45 th Weather Squadron Lightning 
Spreadsheet that Calculates the Probability of any Nearby Lightning Stroke 
being Inside any Radius of Any Point of Interest 


Conversions Module 

Const pie As Double = 3.14159265358979 

Const ecc As Double = 0.081819190842622 'eccentricity (e A ) 

Const kmPerMile As Double = 1.852 


Function toRad(Degrees As Double) As Double 
'Converts degrees to radians 
toRad = Degrees * pie / 1 80 
End Function 


Function toDeg(Radians As Double) As Double 
'Converts radians to degrees 
toDeg = Radians *180/ pie 
End Function 


Function applyConfidence(length As Double, newSigma As Double, currentSigma As Double) As 
Double 

'convert length based on confidence interval 
applyConfidence = length * newSigma / currentSigma 
End Function 

Function revertConfidence(length As Double, newSigma As Double, currentSigma As Double) As 
Double 

'convert length based on confidence interval 
revertConfidence = length * currentSigma / newSigma 
End Function 


Function kmToNmi(kilometers As Double) As Double 
'Converts kilometers to nautical miles 
kmToNmi = kilometers / kmPerMile 
End Function 


Function geodToECEF(Lat As Double, Lon As Double) As Variant 
'convert lat-lon to EFG 
Dim ret(2) As Double 
ret(0) = getE(Lat, Lon) 
ret(l) = getF(Lat, Lon) 
ret(2) = getG(Lat) 
geodToECEF = ret 
End Function 
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Function EcefToGeod(E As Double, F As Double, g As Double) As Variant 
'convert EFG coordinate to lat-lon \ 

Dim ret(l) As Double 'array of (lat,lon) 
ret(O) = Atn(g / ((1 - ecc A 2) * Sqr(E A 2 + F A 2))) *180/ pie 
ret(l) = Atn2(E, F) * 180 /pie 
EcefToGeod = ret 
End Function 


Functions Module 

Const pie As Double = 3.14159265358979 

Const eRad As Double = 3443.920086 ’radius of earth in nautical miles 
Const ecc As Double = 0.081819190842622 'eccentricity (e A ) 


Public Function Atn2(ByVal X As Double, ByVal Y As Double) As Double 
Dim at As Double 

If X = 0 Then 
If Y < 0 Then 
Atn2 = -pie / 2 
Elself Y > 0 Then 
Atn2 = pie / 2 
Else 

Atn2 = 0 
End If 

Elself Y = 0 Then 
If X > 0 Then 
Atn2 = 0 

Elself X < 0 Then 
Atn2 = pie 
End If 
Else 

at = Atn(Y / X) 

IfX<0 And Y<0 Then 
Atn2 = at - pie 

Elself X < 0 And Y > 0 Then 
Atn2 = at + pie 
Else 

Atn2 = at 
End If 
End If 

End Function 


Function arcCos(X As Double) 
arcCos = Atn((l - X A 2) A 0.5 / X) 
End Function 


Function addPoints(pl() As Double, p2() As Double) As Variant 
'add two EFG coordinates together f 

Dim ret(2) As Double 
ret(O) = pi (0) + p2(0) 
ret(l) = pl(l) + p2(l) 
ret(2) = pi (2) + p2(2) 
addPoints = ret 
End Function 


Function getEPointDiffs(eLat As Double, eLon As Double, A As Double, B As Double, phi As Double, t 
As Double) As Variant 

'get EFG differentials for a point on ellipse at t radians 

Dim ret(2) As Double 

Dim xl As Double, yl As Double 

xl = getX(A, B, 0, 90 - phi, t) 

yl = getY(A, B, 0, 90 - phi, t) 

ret(0) = -xl * Sin(eLon * pie / 180) - yl * Sin(eLat * pie / 180) * Cos(eLon * pie / 180) 
ret(l) = xl * Cos(eLon * pie / 1 80) - yl * Sin(eLon * pie / 1 80) * Sin(eLat * pie / 1 80) 
ret(2) = yl * Cos(eLat * pie / 180) 
getEPointDiffs = ret 
End Function 


Function getFocalPointDiffs(eLat As Double, eLon As Double, A As Double, B As Double, phi As 
Double) As Variant 

'get EFG differentials for focal point of ellipse 
Dim ret(2) As Double 

Dim fociDist As Double, xl As Double, yl As Double 
fociDist = Sqr(A A 2 - B A 2) 

xl = fociDist * Cos(toRad(90 - phi)) 
yl = fociDist * Sin(toRad(90 - phi)) 

ret(O) = -xl * Sin(eLon * pie / 1 80) - yl * Sin(eLat * pie / 1 80) * Cos(eLon * pie / 1 80) 
ret(l) = xl * Cos(eLon * pie / 180) - yl * Sin(eLon * pie / 180) * Sin(eLat * pie / 180) 
ret(2) = yl * Cos(eLat * pie / 1 80) 

getFocalPointDiffs = ret 
End Function 


Function getDistGeod(latl As Double, lonl As Double, lat2 As Double, lon2 As Double) As Double 
'spherical 

Dim rLatl As Double, rLat2 As Double, rLonl As Double, rLon2 As Double 

rLatl = toRad(latl) 
rLat2 = toRad(lat2) 
rLonl =toRad(lonl) 
rLon2 = toRad(lon2) 
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getDistGeod = getN((latl + lat2) / 2) * arcCos(Cos(rLatl) * Cos(rLonl) * Cos(rLat2) * Cos(rLon2) + 
Cos(rLatl) * Sin(rLonl) * Cos(rLat2) * Sin(rLon2) + Sin(rLatl) * Sin(rLat2)) / 180 * pie 

End Function 


Function greatCircle(latl As Double, lonl As Double, lat2 As Double, lon2 As Double) As Double 
'great circle distance in NM 

Dim rLatl As Double, rLat2 As Double, rLonl As Double, rLon2 As Double 

rLatl = toRad(latl) 
rLat2 = toRad(lat2) 
rLonl =toRad(lonl) 
rLon2 = toRad(lon2) 

greatCircle = (eRad) * arcCos(Cos(rLatl) * Cos(rLat2) * Cos(rLon2 - rLonl) + Sin(rLatl) * 
Sin(rLat2)) 

End Function 


Function Haversine(latl As Double, lonl As Double, lat2 As Double, lon2 As Double) As Double 
'Calculates distance between two latitude and longitude points by the haversine formula in NM 
Dim dlat As Double, dlon As Double, C As Double, A As Double 


dlat = toRad(latl) - toRad(lat2) 
dlon = toRad(lonl) - toRad(lon2) 

A = Sin(dlat / 2) * Sin(dlat / 2) + Cos(toRad(latl)) * Cos(toRad(lat2)) * Sin(dlon / 2) * Sin(dlon / 2) 
C = 2 * Atn2(Sqr(l - A), Sqr(A)) 

Haversine = eRad * C 

End Function 

Function getX(A As Double, B As Double, H As Double, phi As Double, t As Double) As Double 
'get the X cartesian coordinate of a point on the ellipse at t radians 
getX = H + A * Cos(t) * Cos(phi * pie / 1 80) - B * Sin(t) * Sin(phi * pie / 1 80) 

End Function 


Function getY(A As Double, B As Double, k As Double, phi As Double, t As Double) As Double 
'get the Y cartesian coordinate of a point on the ellipse at t radians 
getY = k + B * Sin(t) * Cos(phi * pie / 1 80) + A * Cos(t) * Sin(phi * pie / 1 80) 

End Function 


Function getN(Lat As Double) As Double 

'get radius of the earth as a given latitude, adjusted for eccentricity 
getN = eRad / (1 - ((ecc A 2) * Sin(Lat * pie / 1 80) A 2)) A 0.5 
End Function 


Function getE(Lat As Double, Lon As Double) As Double 
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I 


'get E coordinate from lat-lon 

getE = getN(Lat) * Cos(Lat * pie / 1 80) * Cos(Lon * pie / 1 80) 
End Function 


Function getF(Lat As Double, Lon As Double) As Double 
'get F coordinate from lat-lon 

getF = getN(Lat) * Cos(Lat * pie / 180) * Sin(Lon * pie / 180) 
End Function 


Function getG(Lat As Double) As Double 
'get G coordinate from lat 

getG = getN(Lat) * (1 - ecc A 2) * Sin(Lat * pie / 1 80) 
End Function 


Function getEPoint(ecLat As Double, ecLon As Double, A As Double, B As Double, phi As Double, t As 
Double) As Variant 

'get lat-lon of a point on ellipse at t radians 

Dim ret(l) As Double 

Dim eC() As Double 

Dim point() As Double 

Dim diffs() As Double 

eC = geodToECEF(ecLat, ecLon) 

diffs = getEPointDiffs(ecLat, ecLon, A, B, phi, t) 

point = addPoints(eC, diffs) 

getEPoint = EcefToGeod(point(0), point(l), point(2)) 

End Function 


Function ePointLat(ecLat As Double, ecLon As Double, A As Double, B As Double, phi As Double, t As 
Double) As Double 

'get lat of a point on ellipse at t radians 
Dim geod() As Double 
geod = getEPoint(ecLat, ecLon, A, B, phi, t) 
ePointLat = geod(0) 

End Function 

Function ePointLon(ecLat As Double, ecLon As Double, A As Double, B As Double, phi As Double, t 
As Double) As Double 

'get Ion of a point on ellipse at t radians 
Dim geod() As Double 
geod = getEPoint(ecLat, ecLon, A, B, phi, t) 
ePointLon = geod( 1 ) 

End Function 


Function getAzimuth(latl As Double, lonl As Double, lat2 As Double, lon2 As Double) As Double 
'get azimuth between two points 
Dim angle As Double 
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angle = (Atn2(Cos(toRad(latl)) * Sin(toRad(lat2)) - Sin(toRad(latl)) * Cos(toRad(lat2)) * 
Cos(toRad(lon2) - toRad(lonl)), Sin(toRad(lon2) - toRad(lonl)) * Cos(toRad(lat2))) * 180 / pie) 
If (angle = 0) Then 
getAzimuth = 360 
Else 

getAzimuth = (angle + 360) Mod 360 
End If 

End Function 


Module 2 Code 

Dim point As New Strike 

Const password As String = "PAEc0d3m0nk3y$" 


Sub POI_Change() 
Calculate 
End Sub 


Sub calculateRow() 

Application.EnableCancelKey = xlDisabled 

Worksheets("Results").protect Contents :=False, passWord:=passWord 

Dim ecefO As Double 

Dim interestPoint(l) As Double ' (lat,lon) 

Dim rowNumber_data As Double 
Dim rowNumber_result As Double 
Dim rowData As Variant 
Dim AOI(3) As Double 
Dim lastRow As Double 

If Worksheets("Results").FilterMode = True Then 
Worksheets("Results").ShowAllData 
End If 

, lastRow = Worksheets("Results'').Cells(Rows.Count, "A").End(xlUp).row 
If (lastRow > 6) Then 

Worksheets("Results").Range("A7:R" & lastRow). Select 
Selection.Delete 
End If 

AOI(O) = Range("latLower") 

AOI(l) = Range("latUpper") 

AOI(2) = Range("lonLower") 

AOI(3) = Range("lonUpper") 

interestPoint(O) = Range("poiLat") 

interestPoint(l) = Range("poiLon") 

ecef = geodToECEF(interestPoint(0), interestPoint(l)) 



Call point.setStatics(interestPoint, ecef, Range("newSigma"), Range("currentSigma"), 
Range("radius")) 

Dim xlCalc As XlCalculation 
xlCalc = Application.Calculation 
Application.Calculation = xlCalculationManual 
' On Error GoTo CalcBack 

lastRow = Worksheets("Data").Cells(Rows.Count, "A").End(xlUp).row 

rowNumber_result = 7 

For rowNumber_data = 2 To lastRow 

rowData = Worksheets("Data").Range("A" & rowNumber_data & ":I" & rowNumber_data) 

If (rowData(l, 3) > AOI(O) And rowData(l, 3) < AOI(l) And rowData(l, 4) > AOI(2) And 
rowData(l , 4) < AOI(3)) Then 

Call point.init(rowData(l , 3), rowData(l, 4), rowData(l, 6), rowData(l, 7), rowData(l, 8), True) 

Worksheets("Results").Cells(rowNumber_result, 1) = rowData(l, 1) 
Worksheets(''Results").Cells(rowNumber_result, 2) = rowData(l, 2) 
Worksheets("Results").Cells(rowNumber_result, 3) = point.centerAzimuth 
Worksheets("Results").Cells(rowNumber_result, 4) = point.centerRange 
Worksheets("Results").Cells(rowNumber_result, 5) = point.criticalAzimuth 
Worksheets("Results").Cells(rowNumber_result, 6) = point.criticalRange 
Worksheets("Results").Cells(rowNumber_result, 7) = rowData(l, 5) 
Worksheets("Results").Cells(rowNumber_result, 8) = point.islnside 
Worksheets("Results").Cells(rowNumber_result, 9) = rowData(l, 9) 
Worksheets("Results").Cells(rowNumber_result, 10) = point. Prob_simpson 
Worksheets("Results").Cells(rowNumber_result, 11) = point. rangeDifference 
Worksheets("Results").Cells(rowNumber_result, 12) = point.Lat 
Worksheets("Results").Cells(rowNumber_result, 13) = point.Lon 
Worksheets("Results").Cells(rowNumber_result, 14) = point. A 
Worksheets("Results").Cells(rowNumber_result, 15) = point.B 
Worksheets("Results").Cells(rowNumber_result, 16) = point.phi 
Worksheets("Results").Cells(rowNumber_result, 17) = point. criticalLat 
Worksheets("Results").Cells(rowNumber_result, 18) = point.criticalLon 
'Worksheets("Results").Cells(rowNumber_result, 19) = point.Prob_rician 
'Worksheets("Results").Cells(rowNumber_result, 20) = point. Prob_simpson 

rowNumber_result = rowNumber_result + 1 
End If 

Next rowNumber_data 

lastRow = Worksheets("Results")-Cells(Rows.Count, "A").End(xlUp).row 
Worksheets(”Results").Range("A7:A" & lastRow). NumberFormat = "MM/DD/YYYY" 
Worksheets('’Results").Range("B7:B'’ & lastRow). NumberFormat = "HH:mm:ss.000" 
Worksheets("Results").Range("C7:C" & lastRo w). NumberFormat = "000" 
Worksheets("Results").Range("D7:D" & lastRow). NumberFormat = "0.00" 
Worksheets("Results").Range("E7:E" & lastRow). NumberFormat = "000" 
Worksheets("Results").Range("F7:F" & lastRow). NumberFormat = "0.00" 
Worksheets("Results").Range("J7:J" & lastRow) .NumberFormat = "0.000000%" 
Worksheets("Results").Range("K7:K" & lastRow).NumberFormat = "0.00" 
Worksheets("Results").Range("L7:M" & lastRow) .NumberFormat = "0.0000" 
Worksheets("Results").Range("N7:P" & lastRow) .NumberFormat = "0.00" 
Worksheets("Results").Range("Q7:R" & lastRow) .NumberFormat = "0.0000" 
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f Worksheets("Results").Range("S7:T" & lastRow).NumberFormat = "0.000000% 


If (lastRow > 6) Then 

With Worksheets("Results")*Range("H7:H" & lastRow). Select 

Selection.FormatConditions.Add Type:=xlTextString, String:="Yes", _ 

T extOperator : =xlContains 

Selection.FormatConditions(Selection.FormatConditions. Count). SetFirstPriority 

With Selection.FormatConditions(l).Font 
.Color = -16777024 
.TintAndShade = 0 

End With 

With Selection.FormatConditions( 1 ) .Interior 
.PattemColorlndex = xlAutomatic 
.Color = 5263615 
.TintAndShade = 0 

End With 

End With 

End If 

Application. Calculation = xlCalc 

Worksheets("Results").EnableAutoFilter = True 

Worksheets("Results").protect Contents;=True, passWord:=passWord, AllowFiltering:=True, 
AllowSorting:=True, userInterfaceOnly:=True 

MsgBox ("Data Calculation Complete") 

Exit Sub 


CalcBack: 

Worksheets("Results").protect Contents:=True, passWord:=passWord, AllowFiltering:=True, 
AllowSorting:=True, userInterfaceOnly:=True 
Application.Calculation = xlCalc 

End Sub 


Sub openBrowser() 

Dim rowNumber As Double, 

Dim N As Integer 

Dim ecef() As Double 

Dim interestPoint(l) As Double ' (lat,lon) 

Dim showPath As Boolean 
showPath = ActiveSheet. showPath. value 

interestPoint(O) = Range("poiLat") 
interestPoint(l) = Range("poiLon") 

ecef = geodToECEF(interestPoint(0), interestPoint(l)) 


Call point. setStatics(interestPoint, ecef, Range("newSigma"), Range("currentSigma"), 
Range("radius")) 

rowNumber = Range("nRow") 

N = Range("nPoints") 

Call point.init(Cells(rowNumber, 12), Cells(rowNumber, 13), Cells(rowNumber, 14), 
Cells(rowNumber, 15), Cells(rowNumber, 16), False) 

Call point.setURL(N, showPath) 

’Range("urlString") = point.URL 

Set browser = CreateObject("IntemetExplorer. Application") 

If (Len(point.URL) > 1950) Then 

MsgBox ("Your request contains too many characters. Please select a smaller number of perimeter 
points or turn on off the path attribute.") 

Else 

browser.Navigate (point.URL) 
browser. Visible = True 
End If 

End Sub 


Probability Module 

Const pie As Double = 3.14159265358979 

Function RicianIntegral(U As Double, V As Double, m As Integer) 

'This subroutine computes the Rician Integral Pm when u, v and m are given. 

'It is preferred because both the exponentials are outside the summation loops. 

'This formulation gives essentially the same numerical value as the other one. 

Dim tm As Double, rm As Double, sm As Double, sumTm As Double, sumTmSm As Double 

On Error GoTo overflow 

tm= 1 

rm= 1 

sm = 1 

sumTm = 1 

sumTmSm = 1 

For i = 1 To m 

tm = tm * V / (2 * i) 
rm = rm*U/(2*i) 
sm = sm + rm 
sumTm = sumTm + tm 
sumTmSm = sumTmSm + tm * sm 
Next i 

Ricianlntegral = Exp(-V / 2) * sumTm - Exp(-(U + V) / 2) * sumTmSm 
Exit Function 
overflow: 
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Ricianlntegral = 0 
End Function 


Function numericIntegral(Xprime As Double, Zprime As Double, SigmaXp As Double, SigmaZp As 
Double, Xe As Double, Ra As Double, N As Integer) 

Dim xl As Double, zl As Double, rl As Double, x2 As Double, z2 As Double, r2 As Double 
Dim R As Double, sum As Double, dTheta As Double, thisPhi As Double, Q As Double 

DPhi = 2 * pie / N 

sum = 0 
thisPhi = 0 

xl = (Xprime + Ra) * SigmaZp / SigmaXp 
zl = Zprime 

Fori= 1 ToN 
thisPhi = thisPhi + DPhi 

x2 = (Xprime + Ra * Cos(thisPhi)) * SigmaZp / SigmaXp 
z2 = Zprime + Ra * Sin(thisPhi) 
rl = Sqr(xl * xl + zl * zl) 
r2 = Sqr(x2 * x2 + z2 * z2) 

dTheta = Application. WorksheetFunction.Asin((xl * z2 - x2 * zl) / (rl * r2)) 

R = (rl + r2) / 2 

sum = sum + Exp(-(R / SigmaZp) A 2 / 2) * dTheta 
xl=x2 
zl = z2 
Next i 

Q = -sum / (2 * pie) 

If Xe < Ra Then 

numericlntegral = 1 + Q 
Elself Xe = Ra Then 
numericlntegral = 1 / 2 + Q 
Else 

numericlntegral = Q 
End If 

End Function 


Function simpsonIntegral(Xprime As Double, Zprime As Double, SigmaXprime As Double, 
SigmaZprime As Double, Ra As Double, N As Integer) 

This program computes the collision probability between two objects. 

The integration is performed using Simpson’s Rule. 

Dim Q(1000) As Double, P As Double 

Dim A As Double, B As Double, C As Double, D As Double, E As Double, F As Double, W As 
Double 

Dim H As Double, DH As Double 


pi = 3.14159265358979 


B = Sqr(2) * SigmaXprime 
C = Xprime / B 

D = 1 / (2 * Sqr(2 * pi) * SigmaZprime) 

W = Ra A 2 
DH = Sqr(W) / N 

'Compute Q(I) array from 0 to (N-l) 

For j = 0 To (N - 1) 

H =j * DH 

A = Sqr(W - H A 2) 

zl = A/B -C 
z2 = A / B + C 

ErfZl = ErrorFunction(zl) 

ErfZ2 = ErrorFunction(z2) 

E = (H - Zprime) A 2 / (2 * SigmaZprime A 2) 

F = (H + Zprime) A 2 / (2 * SigmaZprime A 2) 


Q(j) = (Exp(-E) + Exp(-F)) * (ErfZl + ErfZ2) 
Next j 
Q(N) = 0 

'Compute the integral by Simpson's Rule 
sum = 0 

For j = 1 To(N- 1) 

sum = sum + (3 - (-1) A j) * QG) 


Next j 

sum = sum + Q(0) + Q(N) 
P = D * sum * DH / 3 
simpsonlntegral = P 
End Function 


Function ErrorFunction(X) 

'This subroutine computes the Error Function when X is given 
Dim Z, SqrtPi, tm, sum, Max, m, ErfX, ErfcX 
SqrtPi = 1.77245385090552 
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Z = Abs(X) 

If Z <= 4 Then 
tm = Z 
sum = tm 

Max = Int(20 * Z) + 1 
Form= 1 To Max 

tm = -tm * Z * Z * (2 * m - 1) / (m * (2 * m + 1)) 
sum = sum + tm 
Next m 

ErfX = sum * 2 / SqrtPi 
ElseIfZ<= 5.736 Then 
tm = 1 
sum = 1 
For m= 1 To 1 

tm = -tm *(2*m-l)/(2*Z A 2) 
sum = sum + tm 
Next m 

ErfcX = sum * Exp(-Z A 2) / (SqrtPi * Z) 

ErfX = 1 - ErfcX 
Else 
ErfX = 1 
End If 

ErrorFunction = Sgn(X) * ErfX 
End Function 


Strike Class Module 


Const pie As Double = 3.14159265358979 

Const eRad As Double = 3443.920086 'radius of earth in nautical miles 
Const ecc As Double = 0.081 819190842622 'eccentricity (e A ) 


Private p_poi_geo() As Double 
Private p_poi_ecef() As Double 
Private p_newSigma As Double 
Private p_currentSigma As Double 

Private p_strike_geo(l) As Double 
Private p major As Double 
Private p_minor As Double 
Private p_phi As Double 


Tat Ion point of interest [lat,lon] 
'EFG point of interest [E,F,G] 
'user defined confidence level 
'50% confidence 

Tat Ion of strike [lat,lon] 

'major axis in NM 
'minor axis in NM 
'stroke orientation in degrees 


Private 

Private 

Private 

Private 

Private 

Private 

Private 

Private 


p a As Double 'scaled major axis in NM 

p_b As Double 'scaled minor axis in NM 

p strike ecefO As Double 'EFG strike [E,F,G] 
p_criticalT As Double 'angle of closest point in radians. 
p_critical_geo() As Double Tat Ion of closest point [lat,lon] 
p_critical_ecef(2) As Double 'EFG of closest point [E,F,G] 
p_aziCenter As Double 'azimuth to strike center 
p_rangeCenter As Double 'distance to strike center 


Private p_aziCritical As Double 'azimuth to closest point 

Private p_rangeCritical As Double 'distance to closest point 

Private p_rangeDifference As Double 

Private p_islnside As String 

Private p_URL As String 

Private p_focil() As Double 

Private p_foci2() As Double 

Private p_prob_rician As Double 
Private p_prob_numeric As Double 
Private p_prob_simpson As Double 

Private p_radius As Double 'radius around point of interest 


Public Sub setStatics(poi_geo As Variant, poi_ecef As Variant, newSigma As Double, currentSigma As 
Double, radius As Double) 
p_poi_geo = poi_geo 
p_poi_ecef = poi_ecef 

p_newSigma = Sqr(-2 * Log(l - newSigma / 100)) 
p_currentSigma = Sqr(-2 * Log(l - currentSigma / 100)) 
p_radius = radius 
End Sub 


Public Sub init(Lat, Lon, major, minor, phi, adjustAxis As Boolean) 
On Error GoTo initErr 


p_strike_geo(0) = Lat 
p_strike_geo(l) = Lon 
p_phi = phi 
If (adjustAxis) Then 

'input major and minor are in KM, and unsealed 
p_major = kmToNmi(CDbl(major)) 
p_minor = kmToNmi(CDbl(minor)) 


If p_major = 0 Then 
p major = 0.05 
End If 

If p_minor = 0 Then 
p_minor = 0.05 
End If 

p_a = applyConfidence(p_major, p_newSigma, p_currentSigma) 
p_b = applyConfidence(p_minor, p_newSigma, p_currentSigma) 


Else 

'input major and minor are in NM and scaled 
p_a = CDbl(major) 
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p_b = CDbl(minor) 

If p_a = 0 Then 
p_a = 0.1 
End If 

If p_b = 0 Then 
p_b = 0.1 
End If , 

p_major = revertConfidence(p_a, p_newSigma, p_currentSigma) 
p_minor = revertConfidence(p_b, p_newSigma, p_currentSigma) 
End If 


p_strike_ecef = geodToECEF(p_strike_geo(0), p_strike_geo(l)) 

Call setT 

p_critical_geo = getEPoint(p_criticalT) 

p_aziCenter = getAzimuth(p_poi_geo(0), p_poi_geo(l), p_strike_geo(0), p_strike_geo(l)) 
p_aziCritical = getAzimuth(pjpoi_geo(0), p_poi_geo(l), p_critical_geo(0), p_critical_geo(l)) 
p_rangeCenter = Haversine(pjpoi_geo(0), p_poi_geo(l), p_strike_geo(0), p_strike_geo(l)) 
p_rangeCritical = Haversine(p_poi_geo(0), p_poi_geo(l), p_critical_geo(0), p_critical_geo(l)) 
p_rangeDifference = p_rangeCenter - prangeCritical 
Dim diffs() As Double 
Dim point() As Double 

diffs = getFocalPointDiffs(p_strike_geo(0), p_strike_geo(l), p_a, p_b, p_phi) 

point = addPoints(p_strike_ecef, diffs) 

p_focil = EeefToGeod(point(0), point(l), point(2)) 

diffs(O) = -diffs(O) 

diffs(l) = -diffs(l) 

diffs(2) = -diffs(2) 

point = addPoints(p_strike_ecef, diffs) 

p_foci2 = EcefToGeod(point(0), point(l), point(2)) 

Call setProbability 
Exit Sub 
initErr: 
asd = 0 
End Sub 


Public Property Get A() 
A = pa 
End Property 


Public Property Get BO 
B = p_b 
End Property 


Public Property Get criticalT() 
criticalT = p_criticalT 

End Property j 
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Public Property Get criticalLat() 
criticalLat = p_critical_geo(0) 
End Property 


Public Property Get criticalLon() 
criticalLon = p_critical_geo(l) 
End Property 


Public Property Get criticalAzimuth() 
criticalAzimuth = p_aziCritical 
End Property 


Public Property Get centerAzimuth() 
centerAzimuth = p_aziCenter 
End Property 


Public Property Get criticalRange() 
criticalRange = prangeCritical 
End Property 


Public Property Get centerRange() 
centerRange = prangeCenter 
End Property 


Public Property Get rangeDifference() 
rangeDifference = prangeDifference 
End Property 


Public Property Get URL() 
URL = p_URL 
End Property 


Public Property Get Lat() 
Lat = p_strike_geo(0) 
End Property 


Public Property Get Lon() 
Lon = p_strike_geo(l) 
End Property 
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Public Property Get phi() 
phi = p_phi 
End Property 


Public Property Get Prob_rician() 
Prob_rician = p_prob_rician 
End Property 


Public Property Get Prob_numeric() 
Prob_numeric = p_prob_numeric 
End Property 


Public Property Get Prob_simpson() 
Prob_simpson = p_prob_simpson 
End Property 


Public Property Get islnside() 

Dim dl As Double, d2 As Double, t0() As Double, d3 As Double, tl() As Double 
dl = Haversine(p_focil(0), p_focil(l), p_poi_geo(0), p_poi_geo(l)) 
d2 = Haversine(p_foci2(0), p_foci2(l), p_poi_geo(0), p_poi_geo(l)) 
tO = getEPoint(O) 
tl = getEPoint(pie) 

da = Haversine(tl(0), tl(l), t0(0), tO(l)) 

If ((dl + d2) < (da)) Then 
islnside = "Yes" 

Else 

islnside = "No" 

End If 

End Property 


Private Sub setT() 

Dim thisT As Double, minDist As Double, thisDist As Double, negDist As Double, posDist As 
Double, retumT As Double, step As Double 
Dim pointLat As Double, pointLon As Double 
Dim point() As Double 

retumT = -pie 

point = getEPoint(retumT) 

minDist = Haversine(p_poi_geo(0), p_poi_geo(l), point(O), point(l)) 

For thisT = -pie To pie Step pie / 4 
point = getEPoint(thisT) 

thisDist = Haversine(p_poi_geo(0), p_poi_geo(l), point(O), point(l)) 

If (thisDist < minDist) Then 
minDist = thisDist 


retumT = thisT 
End If 
Next thisT 
step = pie / 4 

While step > pie / (2 A 1 6) 
point = getEPoint(retumT - step) 

negDist = Haversine(p_poi_geo(0), p_poi_geo(l), point(O), point(l)) 
point = getEPoint(retumT + step) 

posDist = Haversine(p_poi_geo(0), p_poi_geo(l), point(O), point(l)) 
If (negDist < minDist) Then 
retumT = retumT - step 
minDist = negDist 
Elself posDist < minDist Then 
retumT = retumT + step 
minDist = posDist 
End If 

step = step / 2 
Wend 

p_criticalT = retumT 
End Sub 


Private Function getEPoint(t As Double) As Variant 
Dim ret(l) As Double 
Dim point() As Double 
Dim diffs() As Double 

diffs = getEPointDiffs(p_strike_geo(0), p_strike_geo(l), p a, p_b, p_phi, t) 

point = addPoints(p_strike_ecef, diffs) 

getEPoint = EcefToGeod(point(0), point(l), point(2)) 

End Function 


Private Function getCPoint(t As Double) As Variant 
Dim ret(l) As Double 
Dim point() As Double 
Dim diffs() As Double 

diffs = getEPointDiffs(p_strike_geo(0), p_strike_geo(l), p_radius, p_radius, p_phi, t) 

point = addPoints(p_poi_ecef, diffs) 

getCPoint = EcefToGeod(point(0), point(l), point(2)) 

End Function 


Public Sub setURL(nPoints As Integer, showPath As Boolean) 

Dim thisT As Double, tStep As Double 
Dim point() As Double 

p_URL = "http ://maps. google. com/staticmap?size=640x640&maptype=satellite&markers 
tStep = 2 * pie / nPoints 
For i = 1 To nPoints 
thisT = -pie + i * tStep 
point = getEPoint(thisT) 
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p_URL = p_URL & Round(point(0), 5) & & Round(point(l), 5) & & "tinyyellow|" 

'p_URL = p_URL & point(O) & & point(l) & & "tinyyellow|" 

Next i 

tStep = 2 * pie / 10 
For i = 1 To 10 

thisT = -pie + i * tStep 
point = getCPoint(thisT) 

p_URL = p_URL & Round(point(0), 5) & & Round(point(l), 5) & & "tinywhite|" 

'p_URL = p_URL & point(O) & & point(l) & & "tinywhite|" 

Next i 

p_URL = p_URL & p_critical_geo(0) & & p_critical_geo(l) & & "smallpurple|" 

'p_URL = p_URL & p_foci 1(0) & & p_foci 1(1 )&","& "smallblue|" 

'p_URL = p_URL & p_foci2(0) & & p_foci2(l) & & "smallbluej" 

p_URL = p_URL & p_poi_geo(0) & & p_poi_geo(l) & & "smallgreen|" 

p_URL = p_URL & p_strike_geo(0) & & p_strike_geo(l) & & "smallred" 

If (showPath) Then 

p_URL = p_URL & "&path=" 
tStep = 2 * pie / nPoints 
For i = 1 To nPoints + 1 
thisT = -pie + i * tStep 
point = getEPoint(thisT) 

p_URL = p_URL & Round(point(0), 5) & & Round(point(l), 5) & "|" 

Next i 
End If 

p_URL = p_URL & 

"&sensor=false&key=ABQIAAAAdwjRqRR8FuOdpAOoimTJCBSOxDK051wxOGB6dfDkgLOxwdqZC 

hSForDLWhvadXUn6EEI6WZWYt853w" 

End Sub 


Public Sub setProbability() 

On Error GoTo probErr 

’This program computes the probability of of a lightning strike occurring within a 
’specified distance of an asset of interest. 

’The cross-section is a circle. 

Dim alpha As Double, LonP As Double, LatP As Double, LonS As Double, LatS As Double 
Dim Xprime As Double, Zprime As Double, SigmaXp As Double, SigmaZp As Double, sigma As 
Double 

Dim X As Double, Z As Double, theta As Double 

Dim U As Double, V As Double 

Dim m As Integer, AspectRatio As Double 

p_prob_rician = 0. 
p_prob_numeric = 0 

If (p_major = 0 Or p_minor = 0) Then 
Exit Sub 
End If 

LatP = toRad(p jpoi_geo(0)) ’Latitude of POI (radians) 
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LonP = toRad(p_poi_geo( 1 )) 'Longitude of POI (radians) 

LatS = toRad(p_strike_geo(0)) 'Latitude of Strike Spot (radians) 

LonS = toRad(p_strike_geo(l )) 'Longitude of Strike Spot (radians) 

alpha = toRad(p_phi) 'Heading of Semi -Major Axis of 50% Confidence Ellipse (radians) 


AspectRatio = p_major / p_minor 
X = (LonP - LonS) * Cos(LatS) 

Z = LatP - LatS 


theta = alpha - ((pie / 2) - Atn2(X, Z)) 


Xprime = p_rangeCenter * Cos(theta) 

Zprime = p_rangeCenter * Sin(theta) 

SigmaXp = p_maj or / 1 .17741 0023 

SigmaZp = p_minor / 1 .177410023 

sigma = Sqr(SigmaXp * SigmaZp) 

V 

'Compute Rician Integral 
U = (p_radius / sigma) A 2 

'V = Xprime A 2 / SigmaXp A 2 + Zprime A 2 / SigmaZp A 2 
'm = Int(Sqr(U * V)) + 1 
'p_prob_rician = RicianIntegral(U, V, m) 

'skip numeric method if aspect ratio is low 
'If (m < 25 And AspectRatio < 1 0) Then 
'Exit Sub 
'End If 

'Use numerical method for computing lightning strike probability: 

'p_prob_numeric = numericIntegral(Xprime, Zprime, SigmaXp, SigmaZp, p_rangeCenter, p_radius, 

1000) 

p_prob_simpson = simpsonIntegral(Xprime, Zprime, SigmaXp, SigmaZp, p_radius, 200) 

Exit Sub 
probErr: 

End Sub 
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